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In  earlier  work  [J.  Opt.  Soc.  Am.  A  21, 13-23  (2004)],  we  showed  that  a  combination  of  linear  models  and 
optimum  Gaussian  sensors  obtained  by  an  exhaustive  search  can  recover  daylight  spectra  reliably  from 
broadband  sensor  data.  Thus  our  algorithm  and  sensors  could  be  used  to  design  an  accurate,  relatively 
inexpensive  system  for  spectral  imaging  of  daylight.  Here  we  improve  our  simulation  of  the  multispectral 
system  by  (1)  considering  the  different  kinds  of  noise  inherent  in  electronic  devices  such  as  charge- 
coupled  devices  (CCDs)  or  complementary  metal-oxide  semiconductors  (CMOS)  and  (2)  extending  our 
research  to  a  different  kind  of  natural  illumination,  skylight.  Because  exhaustive  searches  are  expensive 
computationally,  here  we  switch  to  a  simulated  annealing  algorithm  to  define  the  optimum  sensors  for 
recovering  skylight  spectra.  The  annealing  algorithm  requires  us  to  minimize  a  single  cost  function,  and 
so  we  develop  one  that  calculates  both  the  spectral  and  colorimetric  similarity  of  any  pair  of  skylight 
spectra.  We  show  that  the  simulated  annealing  algorithm  yields  results  similar  to  the  exhaustive  search 
but  with  much  less  computational  effort.  Our  technique  lets  us  study  the  properties  of  optimum  sensors 
in  the  presence  of  noise,  one  side  effect  of  which  is  that  adding  more  sensors  may  not  improve  the  spectral 
recovery.  ©  2005  Optical  Society  of  America 

OCIS  codes:  010.1290,  150.0150,  150.2950,  040.0040. 


1.  Introduction 

No  one  nowadays  would  attempt  to  analyze  naked- 
eye  atmospheric  phenomena  such  as  rainbows,  halos, 
glories,  or  coronas  without  using  instruments.  Day¬ 
light  and  skylight  spectra,  for  example,  are  normally 
measured  with  spectroradiometers,  which  are  com¬ 
plex  and  expensive  instruments  that  provide  only  one 
spectrum  per  measurement.  Thus  when  one  mea¬ 
sures  skylight,  the  illumination  arrives  from  either  a 
small  or  large  angular  subtense  of  the  sky,  depending 
on  the  instrument  field  of  view.  Because  researchers 
could  benefit  from  high-resolution  angular  maps  of 
skylight’s  spectral  power  distribution  (SPD),  we  need 
to  measure  many  skylight  spectra  simultaneously 
across  the  sky  dome.  Multispectral  imaging  systems 
let  us  do  so. 

In  recent  years  the  development  and  design  of  mul- 
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tispectral  color-image  acquisition  devices  has  re¬ 
ceived  much  attention  in  color  science.  By  extending 
our  past  research  on  sky  color,1-4  we  offer  here  a 
theoretical  optimum  design  of  a  multispectral  system 
with  3-5  channels  that  can  recover  the  SPD  of  the 
skylight  incident  on  it.  Our  optimum  multispectral 
system  must  estimate  the  spectral  skylight  radiance 
at  each  pixel  of  the  image  based  on  the  response  of  the 
system’s  channels.  This  is  a  classic  inverse  problem 
that  requires  a  mathematical  estimation  algorithm. 
We  opt  for  a  linear  model  based  on  a  principal  com¬ 
ponents  analysis  (PCA)  described  in  Section  2. 

In  the  future,  we  plan  to  build  a  multispectral  sys¬ 
tem  for  imaging  skylight  by  using  liquid-crystal  tun¬ 
able  filters  (LCTFs)  and  a  linear  monochrome  camera 
with  high  dynamic  range.  Therefore  our  simulated 
system  must  model  the  behavior  of  the  optimum  sen¬ 
sors  realistically,  including  the  effects  of  sensor  noise 
(see  Section  3). 

Any  search  for  an  optimum  set  of  Gaussian  sensors 
(those  whose  spectral  sensitivities  are  Gaussian  func¬ 
tions  of  wavelength)  intended  to  recover  skylight 
spectra  with  a  multispectral  imaging  system  depends 
on  several  factors:  the  spectral  response  of  its  sen¬ 
sors,  the  type  and  number  of  sensors,  and  the  noise 
that  always  affects  any  electronic  device.  To  include 
all  these  factors  in  an  exhaustive  search  is  a  highly 
demanding  computational  task.  Our  alternative  ap- 
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proach  greatly  reduces  computing  time  with  a  simu¬ 
lated  annealing  algorithm  that  minimizes  a  cost 
function  (see  Section  4).  To  this  end,  we  propose  a 
single  cost  function  that  evaluates  the  quality  of  our 
recovered  skylight  spectra:  the  colorimetric  and  spec¬ 
tral  combined  metric  or  CSCM  function. 

In  Section  5  we  present  the  results  of  both  recov¬ 
eries  for  two  different  sets  of  skylight  spectra.  One 
dataset  was  measured  in  Granada,  Spain,  and  was 
used  as  a  training  set  both  in  the  PCA  and  in  the 
simulated  annealing  algorithm.  The  second  dataset 
was  measured  in  Owings,  Maryland,  and  was  used  to 
test  the  accuracy  of  our  simulated  spectral  imaging 
system. 

2.  Estimating  Skylight  Spectra  from  Broadband 
Sensor  Data 

The  spectral  response  of  CCD  camera  sensors  is  often 
assumed  to  be  linear.5-7  If  we  make  this  assumption 
for  our  multispectral  imaging  system,  we  can  model 
its  sensors’  response  using 


Pk  —  J  E(\)Rk(\)d\,  (1) 

^min 

where  pk  is  the  response  of  Mh  sensor,  E(\)  is  the 
illuminant  spectrum  (skylight  in  our  case),  and  Rk(\) 
is  the  Mh  sensor’s  spectral  sensitivity.  If  we  sample 
the  visible  spectrum  at  N  different  wavelengths,  Eq. 
(1)  is 

P*  =  2  (2) 

m  =  1 

In  an  earlier  paper,2  we  measured  1567  skylight 
spectra  in  Granada,  Spain  (37°  10'N,  3°  36'  W,  ele¬ 
vation  680  m)  at  many  different  solar  elevations; 
each  spectrum  ranged  from  380  —  780  nm  in  5  nm 
steps.  This  dataset  allowed  us  to  obtain  81  eigenvec¬ 
tors  from  a  principal  components  analysis8-9  (PCA). 
The  eigenvectors  Vn(\m)  can  then  be  used  in  a  linear 
model  to  reconstruct  a  given  skylight  spectrum  by 

E(km)  =  2  (3) 

n  =  l 

This  PCA-based  linear  reconstruction  method  has 
been  used  widely5-6-8-16  and  gives  good  results  com¬ 
pared  with  other  estimation  methods  (such  as  Wie¬ 
ner’s  method,  spline  interpolation,  modified  discrete 
sine  transform,  nonnegative  matrix  factorization,  or 
direct  transformations) . 12- 16- 17 
The  weight  en  of  each  principal  component  in  the 
linear  combination  is  calculated  from  a  camera’s  sen¬ 
sitivity  (or  sensor  response)  spectra.  If  we  substitute 
Eq.  (3)  in  Eq.  (2)  and  express  it  in  matrix  form,  then 


where  p  is  a  vector  of  k  rows  containing  the  k  sensors’ 
responses,  R  is  an  N  X  k  matrix  containing  the  spec¬ 
tral  sensitivities  of  the  k  sensors  at  N  sampled  wave¬ 
lengths  (superscript  T  denotes  its  transpose),  V  is  a 
N  X  n  matrix  containing  the  first  n  eigenvectors  used 
in  reconstructing  N  wavelengths,  and  e  is  a  vector  of 
n  rows  that  contains  the  coefficients  of  Eq.  (3)’s  linear 
combination.  Then  A  is  a  k  X  n  matrix  that  directly 
transforms  the  coefficients  e  into  the  sensor  re¬ 
sponses  p.  By  calculating  A’s  pseudoinverse  (denoted 
by  superscript  +),  we  obtain  the  coefficients  for  the 
linear  estimate  of  the  camera  sensor  responses,  and 
then  we  can  recover  the  skylight  spectrum 

ER(\J  =  VA+p.  (5) 

So  with  a  combination  of  k  sensors  and  n  eigenvec¬ 
tors,  we  can  estimate  skylight  spectra  using  Eq.  (5). 
Note  that  k  and  n  need  not  be  equal,  although  we  find 
that  this  often  gives  the  best  results. 

3.  Influence  of  Noise 

Any  real  imaging  system  is  of  course  affected  by 
noise,5  a  fact  not  explicitly  accounted  for  in  Eqs.  ( 1)— 
(5).  Yet  noise  can  be  represented  there  as  an  additive 
term  that  changes  the  ideal  noise-free  sensor  re¬ 
sponses  to 

Pnoise  P  ~t  13) 

where  or  is  a  Mrow  vector  of  uncorrelated  components 
that  affect  each  sensor  separately.  There  are  various 
sources  of  noise,18  with  thermal  noise  being  the  most 
common.  This  consists  of  random,  normally  distrib¬ 
uted  white  noise  (i.e.,  affecting  every  wavelength 
equally)  that  is  additive  and  whose  variance  in¬ 
creases  with  temperature.  Another  noise  source  in 
electronic  systems  is  shot  noise,  whose  source  is  cur¬ 
rent  fluctuations  in  semiconductor  devices  due  to  the 
quantum  character  of  electrons.  Although  shot  noise 
is  insignificant  compared  with  thermal  noise  in 
metal-oxide  semiconductor  (MOS)  devices  such  as 
CCDs,  it  too  is  normally  distributed.  Flicker  noise  is 
also  common,  and  it  varies  inversely  with  camera 
electric  frequency  (the  temporal  frequency  at  which 
pixels  charges  are  read).  Finally,  quantization  noise 
is  present  in  every  analog-to-digital  (A/D)  conversion 
and  is  the  loss  of  least-significant  digits  when  quan¬ 
tizing  scene  radiances  to  a  given  number  of  bits  (i.e., 
to  a  fixed  number  of  discrete  levels). 

In  this  study,  we  simulate  thermal  and  shot  noise 
as  random,  normally  distributed  noise  with  standard 
deviations  of  1%,  3%,  or  5%  of  the  maximum  sensor 
response  (these  noise  levels  correspond  to  signal-to- 
noise  ratios  (SNRs)  of  40  dB,  30  dB,  and  26  dB,  re¬ 
spectively).  Quantization  noise  is  given  as  that  due  to 
A/D  conversion  at  resolutions  of  8,  10,  or  12  bits.  This 
slightly  degraded  performance  closely  approximates 
the  behavior  of  a  real  multispectral  imaging  system. 


p  =RtVe  =  Ae, 


(4) 
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Table  1.  Means  and  Standard  Deviations  (SD)  for  1567  Skylight  Spectra  Measured  in  Granada,  Spain,  Using  3  Sensors,  3  Eigenvectors,  and  12-Bit 

Quantization  in  Recovering  Spectra  at  Different  Signal-to-Noise  Ratios  (SNR) 


SNR 

GFC  ±  SD 

A E\b  ±  SD 

IIE(%)  ±  SD 

LN(1  +  1000 
(1  -  GFC))  ±  SD 

CSCM  ±  SD 

40  dB 

0.9993  ±  0.0012 

0.7  ±  0.5 

1.3  ±  0.7 

0.4  ±  0.3 

2.4  ±  1.1 

30  dB 

0.9987  ±  0.0016 

0.9  ±  0.5 

3.2  ±  1.9 

0.7  ±  0.4 

4.8  ±  2.2 

26  dB 

0.9981  ±  0.0023 

1.3  ±  0.7 

5.0  ±  3.5 

0.9  ±  0.5 

7.3  ±  4.1 

4.  Searching  for  Optimum  Sensors 

Various  methods  have  been  proposed  for  selecting  the 
optimum  filters  or  sensors  for  a  multispectral  imag¬ 
ing  system. 5“7'14'17’19  However,  no  consensus  exists 
about  what  “optimum”  means  in  such  a  system.  For 
us,  one  set  of  sensors  is  clearly  better  than  another  if 
its  reconstructed  skylight  spectra  are  more  accurate 
when  tested  by  some  metric.  The  key  question  is  what 
metric  to  use. 

For  our  problem,  essentially  two  kinds  of  metrics 
exist:  colorimetric  and  spectral.20  Colorimetric  met¬ 
rics  formulated  by  the  CIE  (Commission  Interna¬ 
tionale  de  l’Eclairage),  include  CIELUV,  CIELAB, 
CIE94,  and  CIEDE2000.21-22  These  metrics  quantify 
distances  in  their  respective  quasi-uniform  color 
spaces  and  approximate  color  differences  observed  by 
the  human  eye.  To  calculate  such  metrics,  we  only 
need  to  know  the  tristimulus  values23  of  a  given  color 
signal.  Note  that  any  skylight  color  can  be  quantified 
by  spectrally  convolving  its  SPD  with  each  of  three 
CIE  color-matching  functions  (we  use  functions  for 
the  2°  1931  CIE  standard  observer).  Integrals  of  the 
three  convolved  spectra  yield  the  particular  skylight 
color’s  tristimulus  values.  This  trichromatic  match¬ 
ing  has  a  major  consequence:  we  cannot  distinguish 
between  metamers,  which  are  color  stimuli  with  the 
same  tristimulus  values  but  different  SPDs.  Spectral 
metrics  are  those  that  measure  the  distance  between 
two  spectral  curves,  such  as  root-mean-square  error 
(RMSE)  or  goodness-fit  coefficient2  (GFC,  which  uses 
Schwartz’s  inequality,  a  widely  accepted20-24  index  of 
similarity  between  two  spectra).  These  metrics  dis¬ 
tinguish  between  metamers  but  do  not  consider  hu¬ 
man  vision.  Some  new  metrics  have  been  proposed  for 
comparing  spectra  that  take  into  account  properties 
of  the  human  visual  system,  such  as  weighted  RMSE 
(WRMSE)  with  the  diagonal  of  Cohen  matrix  R.20 
Viggiano  proposed  a  spectral  comparison  index26 
(SCI),  a  metric  whose  properties  have  been  studied  by 
others.20-24-25  Another  metric  widely  used  in  solar  ra¬ 
diation  measurements  is  the  percentage  of  the  inte¬ 
grated  irradiance  error26  (IIE(%))  across  the  visible 
spectrum. 

Some  authors6  have  searched  for  optimum  sensors 
using  only  one  of  the  metrics  described  above.  Be¬ 
cause  their  results  depend  on  the  metric  used,  they 
are  not  particularly  useful  in  selecting  sensors  for  our 
planned  multispectral  system.  Imai  et  al.  suggest 
that  “mononumerosis”  should  be  avoided  when  eval¬ 
uating  the  quality  of  spectral  matches.20  By  this  term 
they  mean  that  several  metrics  should  be  used  to 


assess  color  reconstruction  from  both  colorimetric 
and  spectral  standpoints.  Day14  used  thresholds  for 
RMSE  and  CIEDE2000  metrics  when  searching  for 
optimum  sensors;  Hernandez-Andres  et  al.1  used 
GFC,  CIELUV,  and  IIE(%)  in  a  similar  way. 

As  explained  below,  we  must  use  a  single  cost  func¬ 
tion  when  developing  a  simulated  annealing  algo¬ 
rithm.  This  approach  may  seem  to  contradict  the 
recommendations  of  Imai  et  al20  Yet  it  does  not,  be¬ 
cause  we  actually  construct  a  simple  single-cost  func¬ 
tion  or  metric  that  combines  several  metrics  at  once. 
We  use  GFC  as  a  spectral  metric,  CIELAB  as  a  col¬ 
orimetric  cost  function  (denoted  by  A E*ab,  the  dis¬ 
tance  between  two  colors  in  the  CIE’s  uniform  color 
space  L*a*b*),  and  IIE(%)  as  a  metric  for  comparing 
the  spectral  curves  of  natural  illuminants.  In  princi¬ 
ple,  our  new  metric  should  approach  zero  for  near¬ 
perfect  matches  (in  contrast  with  GFC,  which  tends 
to  unity  for  perfect  matches)  and  give  approximately 
the  same  weight  to  the  GFC,  CIELAB,  and  IIE(%) 
metrics.  As  defined  in  Section  1,  our  combined  CSCM 
metric  is  calculated  by 

CSCM  =  Ln(l  +  1000(1  -  GFC))  +  A E*ab  +  IIE(%). 

(7) 

Equation  (7)  defines  a  combined  metric  that  is  zero 
for  perfect  matches,  and  thus  is  a  good  candidate  for 
developing  an  annealing  search  algorithm.  The  chief 
advantage  of  this  metric  is  that  it  quantifies  the  fol¬ 
lowing:  spectral  mismatches  among  metamers,  per¬ 
ceptual  differences  in  color  matches,  and  differences 
in  such  integrated  radiometric  quantities  as  irradi¬ 
ance  or  radiance.  Though  this  metric  may  not  avoid 
“mononumerosis,”  it  clearly  combines  the  properties 
of  various  metrics  relevant  to  skylight  spectra. 

Table  1  shows  the  means  and  standard  deviations 
(SDs)  for  our  Granada  skylight  spectra  obtained  us¬ 
ing  the  linear  method  described  above  for  the  3  best 
sensors  (as  described  later)  and  3  eigenvectors  at 
various  noise  levels  (always  with  12-bit  quantiza¬ 
tion).  Note  that  the  GFC,  CIELAB,  and  IIE(%)  terms 
are  roughly  equal  in  each  row  of  Table  1,  thus  justi¬ 
fying  our  weights  for  them  in  Eq.  (7). 

Now  that  Eq.  (7)  quantitatively  defines  our  opti¬ 
mum  sensor,  we  next  turn  to  developing  a  search 
algorithm.  Whenever  possible,  one  should  do  an  ex¬ 
haustive  search  to  find  a  multispectral  system’s  op¬ 
timum  sensors.  Yet  such  a  search  can  demand 
excessive  computer  time  because  the  number  of  pos¬ 
sible  filter  sets  can  be  enormous.  We  perform  our 
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study  using  3-5  sensors  that  are  Gaussian  functions 
of  wavelength.6’7’10"17  These  sensors  are  similar  to 
commercial  ones  and  could  also  be  made  using  a 
LCTF.  We  vary  sensors’  peak  sensitivities  from  380 
-  780  nm  in  5  nm  steps,  a  spectral  resolution  ade¬ 
quate  for  both  colorimetry  and  radiometry.  We  also 
vary  the  sensors’  FWHM  (full  width  at  half  maxi¬ 
mum)  from  10  —  300  nm  in  10  nm  steps.  Finally,  we 
perform  linear  spectral  recoveries  using  3,  4,  and  5 
eigenvectors.  To  appreciate  the  computational  bur¬ 
den  involved,  note  that  ~1010  different  sets  must  each 
be  evaluated  to  find  the  optimum  set  for  a  3-sensor 
system,  a  search  that  requires  several  days  on  exist¬ 
ing  personal  computers.  This  huge  number  grows  if 
we  try  to  find  the  best  4  or  5  sensors,  for  which  the 
task  is  now  unfeasible  because  the  numbers  of  pos¬ 
sible  sets  increase  by  factors  of  103  and  106,  respec¬ 
tively. 

Faced  with  such  daunting  computational  chal¬ 
lenges,  we  turn  to  simulated  annealing  algorithms 
that  greatly  speed  the  finding  of  optimum  solutions  to 
a  system  with  many  different  sets  of  sensors.  If  we 
slightly  relax  our  requirements  for  recovery  accuracy, 
an  annealing  algorithm  can  give  a  nearly  optimum 
solution  after  testing  only  ~105  sets  of  sensors.  In¬ 
deed,  an  annealing  algorithm  gives  ever-better  solu¬ 
tions  the  longer  we  let  it  run,  in  contrast  with  an 
exhaustive  search  in  which  simulation  time  is  not  a 
variable  that  can  be  chosen  a  priori. 

Simulated  annealing  algorithms  have  been  widely 
used  as  search  algorithms  in  physics27-28  and  in  the 
design  of  multispectral  imaging  systems,6’17  but  typ¬ 
ically  they  evaluate  only  a  single  metric  such  as 
CIELAB  A E  ah  or  spectral  RMSE.  Such  algorithms 
are  based  on  simulating  the  process  of  annealing 
(slow  cooling  after  heating)  of  a  thermodynamic  sys¬ 
tem  (e.g.,  a  gas)  that  is  always  in  thermal  equilib¬ 
rium.  The  algorithm  searches  for  the  minimum 
energy  state  when  temperature  decreases  with  time, 
or  at  least  for  a  local  minimum  from  which  the  system 
will  not  move  without  an  enormous  energy  perturba¬ 
tion,  a  condition  not  found  in  thermal  equilibrium. 
We  must  construct  a  rule  for  changing  the  existing 
state,  calculate  its  energy,  and  accept  it  as  the  system 
state  with  a  probability  given  by  Bolztmann’s  factor 
g— ae/at  jn  Qur  cag6;  energy  is  replaced  with  Eq.  (7)’s 
CSCM  cost  function,  and  this  substitution  makes 
clear  why  the  CSCM  must  be  a  single  function  that 
equals  zero  for  perfect  matches.  The  state  is  repre¬ 
sented  by  a  given  set  of  sensors  (the  peak  sensitivities 
and  FWHM  of  which  determine  the  energy  of  the 
state),  and  the  rule  of  state-changing  will  statistically 
favor  those  states  whose  energy  is  close  to  that  of  the 
current  state. 

We  have  compared  the  efficiency  of  the  simulated 
annealing  algorithm  with  an  exhaustive,  yet  feasible, 
search  across  3  sensors.  We  performed  the  exhaustive 
search  together  with  the  simulated  annealing  search 
using  (1)  the  CSCM  as  a  single  cost  function  and  (2) 
various  random  additive  noise  levels,  always  with 
12  bits  for  quantization  noise  (Table  2).  In  most  cases 


the  annealing  algorithm  found  the  same  optimum 
solution  as  the  exhaustive  search.  In  all  other  cases, 
the  two  algorithms’  solutions  were  quite  similar,  and 
the  shapes  of  the  sensors  sensitivity  curves  were 
nearly  identical.  This  suggests  that  the  local  mini¬ 
mum  given  by  the  annealing  algorithm  is  usually 
equivalent  to  that  obtained  by  the  exhaustive  search. 

5.  Results 

In  this  section  we  present  the  best  sensors  obtained 
with  the  simulated  annealing  algorithm  in  various 
cases.  We  first  examine  the  influence  of  the  cost  func¬ 
tion  by  restricting  ourselves  to  just  3  sensors.  Then 
we  use  CSCM  exclusively  to  study  the  influences  of 
noise,  number  of  sensors,  and  number  of  eigenvectors 
on  sensor  spectral  sensitivities.  Our  larger  goal  is  to 
determine  general  properties  of  the  sensors  for  our 
planned  multispectral  imaging  system. 

As  noted  above,  some  authors  use  only  one  metric 
in  their  simulated  annealing  algorithms  to  find  the 
best  sensors  for  a  multispectral  system.6’17  Connah  et 
al.  found  that  using  RMSE  as  a  cost  function  pro¬ 
duced  sensors  evenly  distributed  across  the  visible 
spectrum.6  When  they  used  CIELAB,  they  obtained 
sensors  with  sensitivity  spectra  similar  to  those  of 
human  cones. 

As  Table  3  shows,  using  only  a  single  metric  pro¬ 
duces  sensors  that  work  well  according  to  that  metric 
but  that  perform  poorly  according  to  the  other  met¬ 
rics.  In  particular,  CIELAB  alone  should  not  be  used 
as  a  cost  function  because  its  small  AE*a6  errors  come 
at  the  price  of  large  GFC  and  IIE(%)  errors.  Thus  we 
use  the  CSCM,  which  strikes  a  balance  between  the 
three  different  metrics. 

In  a  second  round  of  simulations,  we  looked  for  the 
best  set  of  3,  4,  and  5  sensors  to  recover  our  Granada 
skylight  spectra  using  the  linear  model  with  3,  4,  and 
5  eigenvectors  at  various  noise  and  quantization  lev¬ 
els. 

Figure  1  shows  that  for  a  given  number  of  sensors, 
the  peak  sensitivities  and  FWHM  are  similar  at  dif¬ 
ferent  noise  levels  for  the  3-5  eigenvectors  used  in  each 
linear  reconstruction  of  a  skylight  spectrum.  This  be¬ 
havior  is  desirable  in  a  practical  multispectral  system. 
As  other  authors  have  noted,7’10  sensor  sensitivity 
curves  tend  to  sharpen  when  the  noise  is  high  (i.e.,  low 
SNR).  This  occurs  because  sharper  sensors  make  the 
matrix  A  more  robust  to  noise  by  decreasing  its  condi¬ 
tion  number.  Not  surprisingly,  the  curves  also  sharpen 
as  we  increase  their  number  (i.e.,  as  we  approach  a 
narrow-band  hyperspectral  imaging  system). 

Table  4  compares  the  mean  values  of  the  metrics 
used  in  this  study  for  our  Granada  skylight  spectra 
recovered  with  the  best  sets  of  sensors.  Both  Table  4 
and  its  graphical  representation  in  Fig.  2  show  that 
for  low  noise  (i.e.,  high  SNR),  recovered  skylight  spec¬ 
tra  are  more  accurate  if  we  ( 1)  increase  the  number  of 
sensors  and  (2)  match  the  number  of  eigenvectors  and 
sensors.15  The  former  effect  occurs  because  we  can 
sample  the  visible  spectrum  more  reliably  with  5  sen¬ 
sors  than  with  3,  although  differences  in  the  metrics 
are  not  large. 
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Table  2.  Comparison  of  Exhaustive  and  Simulated  Annealing  Searches  for  the  Granada  Skylight  Spectra  Recovered  with  3  Sensors  at  12-Bit 

Quantization  for  Different  SNR“ 


SNR 

Search 

Number  of 
Eigenvectors 

Peak  Sensitivities  (nm) 

1st  2nd  3rd 

Sensor  Sensor  Sensor 

1st 

Sensor 

FWHM  (nm) 

2nd 

Sensor 

3rd 

Sensor 

CSCM  ±  SD 

40  dB 

Exhaustive 

3 

380 

465 

615 

250 

80 

210 

2.4  ±  1.1 

Annealing 

380 

460 

630 

180 

80 

120 

2.4  ±  1.1 

Exhaustive 

4 

380 

465 

635 

280 

70 

190 

2.4  ±  1.2 

Annealing 

380 

465 

635 

280 

70 

190 

2.4  ±  1.2 

Exhaustive 

5 

395 

475 

630 

260 

60 

250 

2.5  ±  1.2 

Annealing 

395 

475 

630 

260 

60 

250 

2.5  ±  1.2 

30  dB 

Exhaustive 

3 

380 

465 

595 

170 

90 

120 

4.8  ±  2.2 

Annealing 

380 

465 

595 

170 

90 

120 

4.8  ±  2.2 

Exhaustive 

4 

395 

470 

620 

250 

80 

190 

5.0  ±  2.3 

Annealing 

395 

470 

620 

230 

70 

180 

5.0  ±  2.3 

Exhaustive 

5 

395 

465 

615 

250 

80 

190 

4.9  ±  2.1 

Annealing 

395 

465 

615 

250 

80 

190 

4.9  ±  2.1 

26  dB 

Exhaustive 

3 

395 

460 

550 

100 

60 

90 

7.3  ±4.1 

Annealing 

395 

460 

550 

100 

60 

90 

7.3  ±4.1 

Exhaustive 

4 

385 

460 

605 

200 

70 

150 

7.1  ±  3.2 

Annealing 

385 

460 

605 

200 

70 

150 

7.1  ±  3.2 

Exhaustive 

5 

395 

470 

605 

220 

70 

160 

6.7  ±  3.0 

Annealing 

400 

465 

605 

200 

70 

160 

6.9  ±  2.8 

“Cases  where  the  annealing  and  exhaustive  searches  found  the  same  solution  are  in  bold  type. 


Something  similar  occurs  when  we  quantize  at  ei¬ 
ther  12,  10,  or  8  bits,  as  seen  in  Fig.  3  for  the  case  of 
5  sensors,  5  eigenvectors,  and  an  SNR  of  30  dB.  Re¬ 
sults  improve  if  we  use  12  bits  rather  than  8  bits,  but 
the  difference  is  small.  Thus  we  can  use  a  cheaper 
and  faster  8-bit  A/D  converter  without  significantly 
degrading  system  accuracy.  In  fact,  the  most  impor¬ 
tant  step  in  building  a  multispectral  imaging  system 
is  to  optimize  the  sensors  for  the  remote-sensing  task 
at  hand  rather  than  to  increase  their  number  or 
quantization  levels,  both  of  which  will  increase  sys¬ 
tem  cost  and  response  time. 

For  high-noise  (low  SNR)  sensors,  increasing  the 
number  of  eigenvectors  from  3  to  5  is  always  pref¬ 


erable,  whereas  adding  more  sensors  simply  in¬ 
creases  system  noise.  This  can  be  appreciated  by 
noting  that  a  system  with  more  sensors  likely  has 
more  connections,  more  transistors  in  the  CCD,  and 
more  memory  cells.  All  these  elements  individually 
contribute  noise,  and  so  each  raises  the  total  noise 
level.  As  other  authors  have  noted,  1’7>10  spectral  re¬ 
coveries  do  not  improve  significantly  in  noise-free 
simulations  if  the  number  of  sensors  increases 
from,  say,  four  to  seven  (the  particular  numbers  of 
sensors  depend  on  system  hardware  and  on  the 
shapes  of  the  skylight  spectra).  If  we  add  noise  to 
such  simulations,  using  more  sensors  increases  the 
noise  in  matrix  A.  This  noise  propagates  through- 


Table  3.  Comparison  of  Best  3  Sensors  Found  Using  Annealing  Searches  with  Various  Metrics,  3  Eigenvectors,  and  the  Granada  Skylight  Spectra" 


Peak  Sensitivities  (nm)  FWHM  (nm) 


Cost 

Function 

1st 

Sensor 

2nd 

Sensor 

3rd 

Sensor 

1st 

Sensor 

2nd 

Sensor 

3rd 

Sensor 

GFC 

±  SD 

A E*ab 

±  SD 

IIE(%) 

±  SD 

CSCM 
±  SD 

GFC 

400 

470 

645 

50 

40 

60 

0.9994 

±  0.0012 

0.8  ±  0.6 

1.3  ±  0.7 

2.5  ±  1.1 

A E*ab 

445 

525 

605 

30 

40 

50 

0.9923 

±  0.0131 

0.3  ±  0.2 

2.3  ±  1.6 

4.3  ±  2.3 

IIE(%) 

470 

475 

720 

20 

250 

290 

0.9972 

±  0.0040 

1.3  ±  0.8 

0.7  ±  0.5 

3.2  ±  1.4 

CSCM 

380 

460 

630 

180 

80 

120 

0.9993 

±  0.0012 

0.7  ±  0.5 

1.3  ±  0.7 

2.4  ±  1.1 

“SNR  =  40  dB,  and  quantization  is  at  12  bits.  The  best  result  for  each  metric  appears  in  bold  type,  when  it  alone  was  the  annealing 
algorithm’s  cost  function. 
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wavelength  (nm) 
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-3  sensors 
-4  sensors 
-5  sensors 


Fig.  2.  Means  for  the  CSCM  metric  for  1567  skylight  spectra 
measured  in  Granada,  Spain,  using  3-5  sensors  and  3-5  eigenvec¬ 
tors  at  12-bit  quantization  for  various  SNR. 


wavelength  (nm) 

(C) 


Fig.  1.  Each  row  shows  the  spectral  sensitivity  curves  of  the  best 
3-,  4-,  or  5-sensor  systems,  respectively  (the  numbers  of  eigenvec¬ 
tors  and  sensors  are  equal  in  each  case).  Solid  curves  represent 
sensors  for  SNR  =  40  dB,  while  dashed  curves  represent  sensors 
for  SNR  =  26  dB. 


out  the  mathematical  path  described  in  Section  2 
and  degrades  system  performance  more  than  the 
additional  sensors  improve  it.7-10 

Several  steps  are  needed  to  implement  a  working 
multispectral  system  such  as  those  simulated  here. 
First,  we  must  establish  a  priori  the  accuracy  re¬ 
quired  of  the  system.  Second,  we  must  reduce  noise 
by  all  practical  means  (e.g.,  cooling  the  CCD  and 
subtracting  its  dark  noise).  Finally,  we  must  decide 
how  many  sensors  our  system  will  have,  select  a  pre¬ 
ferred  A/D  converter,  and  calculate  both  the  system’s 
desired  response  times  and  cost. 

Naturally  our  system  must  work  for  spectra  other 
than  those  used  to  calculate  the  original  dataset’s 
eigenvector  matrix  V.  Thus  we  extend  our  analysis  to 
skylight  spectra  measured  at  another  site  in  Owings, 
Maryland.  Using  12-bit  quantization  and  the  best  set 
of  sensors  found  for  each  case  (see  Table  4),  we  ana¬ 
lyze  metrics  for  these  new  spectra  in  bold-type  rows 
in  Table  4.  It  too  shows  small  errors  for  the  recovered 
spectra,  demonstrating  that  both  our  spectral  recov¬ 
ery  method  and  optimum  sensors  can  be  used  to  de¬ 
velop  a  reliable  system  for  imaging  skylight  spectra. 
Note  that  although  mean  GFC  and  CIELAB  A E*ab 
are  similar  for  the  Owings  spectra,  mean  IIE(%)  is 


SNR  30dB 


CSCM 


-■—CIELAB 
“A—  IIE(%) 

— • —  LN(1-1000(1-GFC)) 


Fig.  3.  Effects  of  quantization  noise:  mean  values  for  metrics  of 
the  Granada  skylight  spectra  using  5  sensors  and  5  eigenvectors 
for  SNR  =  30  dB. 
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Table  4.  Means  for  the  Granada  Skylight  Spectral  Recoveries  Using  the  Best  Sets  of  3  to  5  Sensors  for  Various  SNR  Levels  and  12-Bit  Quantization*1 


SNR 

Number  of 
Sensors 

Number  of 
Eigenvectors 

GFC  ±  SD 

A E*ab  ±  SD 

IIE(%)  ±  SD 

CSCM  ±  SD 

40  dB 

3 

3 

0.9993  ±  0.0012 

0.7  ±  0.5 

1.3  ±  0.7 

2.4  ±  1.1 

3 

5 

0.9988  ±  0.0023 

1.0  ±  0.6 

0.9  ±  0.6 

2.5  ±  1.2 

4 

3 

0.9994  ±  0.0011 

0.8  ±  0.6 

1.2  ±  0.7 

2.4  ±  1.2 

4 

4 

0.9997  ±  0.0003 

0.6  ±  0.3 

1.2  ±  0.5 

2.0  ±  0.7 

5 

3 

0.9993  ±  0.0012 

0.8  ±  0.5 

1.1  ±  0.7 

2.3  ±  1.2 

5 

5 

0.9998  ±  0.0002 

0.3  ±  0.2 

1.1  ±  0.8 

1.6  ±  0.8 

5 

5 

0.9990  ±  0.0002 

0.2  ±  0.1 

0.5  ±  0.4 

1.5  ±  0.5 

30  dB 

3 

3 

0.9987  ±  0.0016 

0.9  ±  0.5 

3.2  ±  1.9 

4.8  ±  2.2 

3 

5 

0.9939  ±  0.0067 

1.6  ±  0.7 

1.7  ±  1.4 

4.9  ±  2.1 

4 

3 

0.9990  ±  0.0017 

0.9  ±  0.4 

3.5  ±  2.1 

4.9  ±  2.4 

4 

4 

0.9991  ±  0.0008 

0.9  ±  0.5 

3.5  ±  1.4 

4.9  ±  1.8 

5 

3 

0.9985  ±  0.0017 

1.1  ±  0.6 

2.8  ±  1.6 

4.7  ±  2.0 

5 

5 

0.9991  ±  0.0009 

0.8  ±  0.4 

3.3  ±  2.1 

4.6  ±  2.3 

5 

5 

0.9984  ±  0.0009 

0.7  ±  0.3 

1.7  ±  1.3 

3.3  ±  1.6 

26  dB 

3 

3 

0.9981  ±  0.0023 

1.3  ±  0.7 

5.0  ±  3.5 

7.3  ±  4.1 

3 

5 

0.9867  ±  0.0132 

2.0  ±  0.7 

2.4  ±  2.1 

6.9  ±  2.8 

4 

3 

0.9986  ±  0.0017 

1.2  ±  0.6 

5.9  ±  2.9 

7.9  ±  3.2 

4 

5 

0.9902  ±  0.0014 

2.0  ±  0.8 

2.7  ±  2.2 

6.9  ±  3.1 

5 

3 

0.9976  ±  0.0025 

1.4  ±  0.7 

5.7  ±  3.1 

8.2  ±  3.5 

5 

5 

0.9988  ±  0.0012 

1.1  ±  0.5 

6.1  ±  2.6 

7.9  ±  3.1 

3 

5 

0.9970  ±  0.0018 

1.6  ±  0.7 

1.9  ±  1.4 

4.9  ±  2.6 

“Rows  in  bold  type  are  the  means  for  242  skylight  spectra  measured  in  Owings  recovered  with  the  best  set  of  sensors  at  12-bit 
quantization. 


always  less  for  the  smaller  Owings  dataset  than  for 
the  larger  one  from  Granada.  This  occurs  because  the 
Owings  spectra  are  from  a  limited  range  of  solar  el¬ 
evations  and  are  spectrally  similar  to  most  of  the 
Granada  data.  Both  factors  make  the  Granada  eig¬ 
envectors  well  matched  to  those  that  we  could  calcu¬ 
late  independently  for  Owings.  The  few  spectra 
measured  at  low  solar  elevations  in  Granada  exert 
relatively  little  influence  on  the  PCA  eigenvectors, 
and  these  Granada  spectra  are  the  ones  that  increase 
the  mean  Granada  IIE(%).  Yet  even  the  least  accu¬ 
rate  recoveries  of  skylight  spectra  at  Owings  (see  Fig. 
4)  show  relatively  small  differences  between  the  orig¬ 
inal  and  recovered  spectra. 

6.  Conclusions 

We  have  shown  that  linear  methods  based  on  PCA 
allow  accurate  recovery  of  skylight  spectra  from 
broadband  camera  sensors.  We  propose  CSCM  as  a 
single  metric  that  takes  into  account  three  different 
accuracy  standards:  spectral,  colorimetric,  and  total 
integrated  irradiance.  Our  CSCM  metric  can  easily 
be  used  instead  of  CIELAB  or  RMSE  alone  in  search 
algorithms.  We  have  presented  a  simulated  anneal¬ 
ing  algorithm,  using  CSCM  as  a  single  cost  function, 
as  a  fast  method  for  searching  a  limited  number  of 
Gaussian  sensors  to  construct  an  optimum  multi- 
spectral  imaging  system.  We  have  simulated  some 
common  noise  sources  present  in  any  digital  imaging 
system  in  order  to  mimic  noise  in  real  images.  We 
have  shown  that  increasing  the  number  of  sensors 
does  not  necessarily  improve  the  accuracy  of  recov¬ 


ered  spectra  if  sensor  noise  is  high  because  each  sen¬ 
sor’s  individual  noise  contributions  degrade  the 
overall  quality  of  the  spectral  reconstruction.  Thus 
the  optimum  number  of  sensors  depends  on  noise 
levels  inherent  to  the  given  multispectral  system. 

Skylight  has  complicated  spectra  with  different  ab¬ 
sorption  bands  that  depend  on  species  such  as  water 
vapor,  molecular  oxygen,  ozone,  and  aerosols,  and  the 
relative  strength  of  these  bands  varies  daily  and  even 
hourly.  Yet  despite  this  ever-changing  spectral  detail, 
we  find  that  a  linear  recovery  algorithm  using  only  a 


wavelength  (nm) 

Fig.  4.  Original  skylight  spectrum  (solid  curve)  measured  in 
Owings,  Maryland,  and  recovered  (dotted  curve)  spectrum  for 
SNR  =  40  dB  with  5  sensors,  5  eigenvectors,  and  12-bit  quantiza¬ 
tion  for  the  90th  percentile  of  the  CSCM. 
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few  optimum  Gaussian  sensors  returns  high-quality 
reconstructions  of  skylight  spectra,  even  with  sensor 
noise.  In  future  work,  we  will  use  a  LCTF  and  a 
monochrome  camera  to  build  the  sensors  simulated 
here  in  order  to  study  their  actual  performance  in 
recovering  skylight  spectra. 
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